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A  computer  code,  FEASBL,  is  developed  to  maximize  a 
non-linear  objective  function  over  a  convex  feasible  region. 
The  feasible  region  is  defined  by  a  set  of  non-linear  and 
linear  constraints  on  the  variables.   FEASBL  can  solve  prob- 
lems involving  up  to  fifty  variables  with  a  feasible  region 
formed  by  up  to  fifty  non-linear  constraints,  and  fifty 
linear  constraints.   FEASBL  uses  a  feasible  direction  method 
as  its  solution  algorithm. 
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I.   INTRODUCTION 

This  thesis  presents  a  computer  code  for  solving  medium 
sized  non-linear  programming  problems,  and  describes  its  use 
at  the  Naval  Postgraduate  School  computer  center. 

The  general  non-linear  programming  problem  is  the  optimi- 
zation of  an  objective  function  over  a  constrained  feasible 
region.   By  the  technique  of  separable  programming,  many  non- 
linear problems  can  be  converted  to  a  linear  problem  and 
solved  using  established  linear  programming  methods.   Many 
problems,  however,  cannot  be  solved  by  this  method.   Several 
techniques  have  been  developed  for  solving  non-linear  prob- 
lems.  Most  methods  involve  following  the  gradient  of  the 
objective  function  to  generate  a  sequence  of  points  that 
converges  to  the  optimal  point.   The  method  used  in  this 
code  is  Zoutendijk's  method  of  feasible  directions  [1].   The 
method  of  feasible  directions  can  be  considered  a  large-step 
gradient  method  as  opposed  to  the  small-step  gradient  and 
search  methods.   Feasible  direction  methods  can  be  charac- 
terized by: 

1.  Starting  at  an  initial  feasible  point. 

2.  Establishing  a  feasible  direction  that  will  improve 
the  objective  function. 

3.  Proceeding  to  a  new  feasible  point  that  has  a  better 
value  of  the  objective  function  than  the  previous  point. 

4.  Converging  to  the  optimal  point,  or  determining  that 
the  problem  is  unbounded. 


Many  procedures  can  be  considered  feasible  directions  methods, 
for  example,  the  simplex  method  for  solving  linear  program- 
ming problems. 

The  method  used  in  this  code  requires  that  the  user  know 
the  objective  function,  constraint  functions,  all  their  first 
partial  derivatives,  and  an  initial  feasible  point.   If  the  user 
does  not  know  an  initial ' feasible  point ,  this'  code  can  be 
used  to  solve  a  preliminary  programming  problem  which  will 
provide  a  feasible  point.   This  problem  will  be  described 
later.   This  code  will  solve  maximization  problems  of  the 


form: 


Max:       z  =  F(x) 


Subj  to:   G-(x)  <  b.  ,   ielr 
1  —     1       >~. 


T 
a-  x   <  b.  ,   ieIT 
1  —  -   1  '     L 


d.  <x.  <c.  ,   jeJ  . 


The  functions  F (x) ,  and  G-(x)  are  non-linear  functions  which 
are  supplied  by  the  user  as  Fortran  FUNCTION  routines.   The 
sets  Ir  and  IT  are  the  sets  of  subscripts  for  the  non-linear, 
and  linear  constraints  respectively.   If  the  objective  func- 
tion is  concave,  and  the  non-linear  constraints  are  convex 
functions,  then  the  problem  will  be  a  convex  programming 
problem,  and  the  method  will  converge  to  the  global  maximum. 
If  the  problem  is  not  convex,  then  the  method  may  converge 
to  a  local  maximum.   The  code  is  designed  to  accommodate 
functions  with  up  to  fifty  variables.   The  constraint  set 
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can  include  up  to  fifty  non-linear  constraints  and  fifty 
linear  constraints.   A  complete  description  of  how  to  use 
this  code  is  provided  in  Section  IV. 


II.   MATHEMATICAL  BACKGROUND   . 

The  theoretical  foundation  for  the  method  of  feasible 
directions  was  developed  by  Zoutendijk  in  1960  [1].   A  brief 
discussion  of  the  mathematical  background  as  it  applies  to 
the  method  used  in  this  computer  code  is  presented. 

A.   DEFINITIONS 

A  region  R  is  a  convex  region  if  for  x-i  and  x~ ,  two  points 
in  R,  then  the  point  x.  =  ax.  +  (1  -  a)x?  is  in  R  for  all  a 
such  that  0  <  a  <  1. 

Let  R  be  the  feasible  region,  and  if  x  is  a  member  of  R 
then  a  direction  s_  is  a  feasible  direction  if  and  only  if 
there  exists  a  X'  >  0  such  that  x  +  As  is  in  R  for  all 
Xe(0,A') . 

A  usable  feasible  direction  is  a  feasible  direction  s_ 

T 
such  that  V  F(x)   s  >  0.   In  other  words,  a  usable  feasible 

direction  at  x  is  a  direction  in  which  the  value  of  the  ob- 
jective function  increases  and  along  which.it  is  possible  to 
move  while  remaining  in  R. 

Let  F(x)  be  a  concave  function,  and  let  the  feasible 
region  R  be  convex.   Then  if  there  exists  no  usable  feasible 
direction  at  x  in  R,  then  x  must  be  a  maximum  of  F(x)  over 
the  region  R.   To  show  this  assume  that  x  is  not  a  maximum 
and  there  is  no  usable  feasible  direction  at  x.   There  exists 
a  point  y_  in  R  such  that  F(y_)  >  F(x)  and  since  R  is  convex 
all  the  points  along  the  line  between  y_  and  x  must  be  in  R. 
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T 
Further,  since  F(  )  is  concave,  V  F(x)  (y_  -  x)  >  0,  thus 

the  direction  v_  -  x  is  a  usable  feasible  direction.   This 

contradiction  implies  that  x  must  be  a  maximum  of  F (x ) . 

B.   THE  METHOD 

The  method  of  feasible  directions  begins  at  an  initial 
feasible  point  x  ,  determines  the  "best"  usable  feasible 
direction,  moves  along  that  direction  to  the  point  x   that 
maximizes  F(x)  along  the  direction  without  leaving  the  fea- 
sible region,  repeats  this  procedure  until  a  point  x*  is 
found  which  has  no  usable  direction,  or  determines  that  the 
problem  is  unbounded.   The  method  can  thus  be  viewed  as  a 
four  step  procedure. 

1.  Determine  a  usable  feasible  direction  s_  at  the 
point  x  . 

2.  If  no  usable  feasible  direction  exists  stop,  x   is 

the  maximum  of  F (x)  over  R. 

k+1    k      k 

3.  Determine  a  step  length  A,  such  that  x    =  x   +  A,s_ 

k+1  k 

is  in  R  and  x    maximizes  F(x)  along  s_  . 

4.  If  A,  goes  to  infinity,  then  the  problem  is  unbounded, 
otherwise  go  back  to  step  1  and  repeat. 

The  direction  finding  problem  will  be  discussed  first, 
and  then  the  step  length  problem. 

1 .   The  Direction  Finding  Problem 

The  direction  finding  problem  is  to  find  the  usable 
feasible  direction  that  will  result  in  the  greatest  apparent 
increase  in  the  objective  function.   This  can  be  formulated 
as  the  following  linear  programming  problem. 
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Max:      a 

Subj  to:   a  -  V  F(x)  s  <  0 

a  +  V^Cx^s  <  0,  i£lc(x) 

T 
a.  s_  <  0 ,  iel-,  (x) 

s.  >  0,  jeJ"(x) 
s.    <   0,  jeJ+(x) 

\s.  |  <  1,  jeJ(x), 

where  I  (x)  and  I, (x)  are  the  sets  of  binding  non-linear  con- 
straints and  binding  linear  constraints,  and  J  (x)  and  J  (x) 
are  the  sets  of  j  such  that  x.  is  at  its  lower  and  upper 
bound  respectively. 

If  there  exists  at  least  one  point  that  is  strictly 
interior  to  R  with  respect  to  the  non-linear  constraints  and 
if  the  optimal  value  of  a  is  0  in  the  direction  finding  prob- 
lem, then  there  is  no  usable  feasible  direction  and  we  are 
at  the  optimum  point  for  the  non-linear  problem.   The  last 
constraint  in  the  direction  finding  problem  is  a  normalization 
constraint.   In  order  to  insure  that  the  procedure  will  con- 
verge  the  components  of  s   must  be  bounded,  and  several  nor- 
malizations have  been  studied  by  Zoutendijk.   These  normalizations 
are : 

Nl:    s_Ts  <  1 
N2:   Is. I  <  1 
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N3:    s.  <  1   if  V  F(x)  .  >  0 
J  x    J 

s.  >-l    if  V  F(x) .  <  0 
J  -         x  ~  J 

N4:    a   <    1 

N5:    VxGi(x)T  s  +  a   <   b.  -  G^x);   iel^x) 

T  T 

a-  s  <  b .  -a.  x;  iel. (x) 
— i  —  -   l  — i  — '     L  — J 

d.  -  x.  <  s .  <  c .  -x.. 
J     J  "   3  "   J     J 

The  normalization  Nl  will  generally  lead  to  fewer  steps  then 
N2 ,  N3,  or  N4 ,  but  the  amount  of  computation  for  each  step 
is  considerably  greater.   Also  the  simple  constraints  s.  <  0, 
and  s.  >  0  increases  the  size  of  the  direction  finding  prob- 
lem if  Nl  is  used,  but  not  for  N2 ,  N3,  and  N4.   As  one  goes 
from  N2,  to  N4  the  number  of  steps  needed  increases,  but  the 
amount  of  computation  for  each  step  decreases.   N5  leads  to 
a  larger  problem,  and  also  does  not  converge  as  fast  as  the 
other  normalizations.   Normalization  N2  was  chosen  for  this 
code  as  the  best  compromise  between   the  number  of  steps  and 
the  amount  of  work  for  each  step. 

2 .   The  Modified  Direction  Finding  Problem 

To  insure  that  the  procedure  actually  converges  to 
the  maximum,  some  precautions  must  be  taken  to  prevent  the 
phenomenon  of  "jamming"  in  which  the  procedure  in  effect 
gets  "stuck  in  a  corner"  and  never  finds  the  real  maximum. 
This  is  accomplished  by  considering  nearby  constraints  in 
addition  to  the  binding  constraints  when  solving  the  direc- 
tion finding  problem. 
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Define  the  sets: 

Ir(x,e)  =  {ielr  :  b.  -  e  <  G. (x)  <  b.} 

T 
IT  (x,e)  =  {ieIT  :  b.  -  e  <  a.  x   <  b.} 
L  — '  *  L    i     -—l—   -l 

j"(x,e)  =  {je:J   :d.      <  x.  <  d.  +  e} 

J  (x,e)  =  (jeJ   :  c .  -  e  <  x .  <  c. } . 

The  direction  finding  problem  is  then  modified  to  include  the 
constraints  in  the  sets  Ir(x,e),  I,(x,e),  J  (x,e),  and  J  (x,£) 
Therefore  the  direction  finding  problem  becomes: 

Max :      o 

Subj  to:  a    -    V  F(x)Ts_  <  0 

-A. 

+  VxGi(^)T  -  °  '  i£lc(?_»£) 

T 
a.  s  <  0  ,  ieIT  (x,e) 

s.   >  0  ,  jeJ  (x,e) 

s.   <  0  ,  jeJ  (x,e) 


Is^l  <  1  ,  jeJ(x) 


Then  if  the  solution  to  the  direction  finding  problem  is  such 
that  a  <  e ,  the  problem  is  solved  again  with  e/2.   Thus  we 
stay  away  from  the  "corners"  until  we  get  near  the  optimal 
point . 

3.   The  Step  Length  Problem 

Once  the  direction  s_  has  been  found,  the  problem 

is  to  determine  how  far  to  move  in  that  direction.   If  A.„ 

k     k 
is  defined  to  be  the  largest  step  length  X  such  that  x   +  As_ 


14 


k      k 
is  in  R,  then  the  step  length  will  be  X,  such  that  F(x   +  X,s_  ) 

k     k 
is  the  maximum  of  F(x   +  Xs_  )  :  X  <   X„.   If  X,     =   °°,  then  the 

problem  is  unbounded,  and  we  can  conclude  that' the  maximum  of 

F(x)  is  infinite.   The  maximum  feasible  step  Xp  is  found  by 

k  k 

a  simple  search  along  the  direction  s_  from  the  point  x  . 

Then  the  optimal  step  length  X,  is  determined  by  finding  the 

k     k   k 
X  such  that  V  F(x   +   s  )  s   =0.   Since  F(x)  is  concave,  if 

k      k  T   k 
V  F(x   +  XFs_  )   s_  >  0 ,  then  Xp  is  used  for  the  step  length, 

otherwise  X,  can  be  found  by  vegula   falsi    [2] . 
4.   Stopping  Conditions 

The  procedure  will  terminate  whenever  the  solution 
to  the  direction  finding  problem  yields  a   =  0,  indicating 
that  there  are  no  usable  feasible  directions.   If  the  solu- 
tion to  the  step  length  problem  is  X,  =  °°,  then  the  problem 
is  unbounded  and  the  procedure  is  terminated  indicating  an 
unbounded  solution.   In  addition,  if  the  step  length  or  the 
percent  improvement  in  the  objective  function  becomes  smaller 
than  user  provided  limits,  the  problem  is  terminated,  and 
the  current  solution  is  specified. 

C.   THE  INITIAL  FEASIBLE  POINT 

If  an  initial  feasible  point  is  not  known,  pick  an  x' 
such  that  d.  <  x!  <  c,  and  solve  the  problem: 

Max:       -z1 


Subj  to:     T  .  ,   .  .  T 

J  a.  x  -  p.  z-  <  d.  :  iel 


-i--   t-i"! 


d .  <  x .  <  c  . 

J    -      J  "   J 
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where 


zl   I   ° 


Pi  = 


0  ;  if  a-  x'  <  b. 
1  —  -   1 


a-  x'  -  b.  ;  if  a.  x'  >  b. , 
— l  —     i      — r  —     i 


An  initial  feasible  point  for  this  linear  programming  problem 
is  x1,  z.,  =  1.   If  the  feasible  region  is  non-empty,  then  the 
optimal  solution  of  this  problem  will  be  z,  =  0,  x  =  x  " .   Now 
x"  will  be  feasible  with  respect  to  the  linear  constraints, 
if  x"  is  also  feasible  with  respect  to  the  non-linear  con- 
straints, then  use  x"  for  the  initial  feasible  point.   How- 
ever if  x"  for  the  initial  feasible  point.   However  if  x"  is 
hot  feasible  with  respect  to  the  non-linear  constraints, 
solve  the  problem: 


Max: 


Subj  to: 


G-  (x)  -  p.  z„  <  b .  ;  ielr 
*x  2  -   l  '     C 

<  b-  ;  ieIT 


T 


d.  <  x.  <  c . 
J  "   J  -   J 

z2  >  0 


where 


Pi  = 


;  if  G. (x")  <  b,  ) 


G. (x")  -  b.  ;  if  G.  (xM)  >  b.  . 
l  —  *  i       l  —  *  l   / 
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This  non-linear  programming  problem  can  be  solved  by  the 
method  of  feasible  directions  using  x  " ,  z-  =  1  as  the  initial 
feasible  point.   If  the  feasible  region  of  the'  original  prob- 
lem is  non-empty,  then  the  optimal  solution  to  this  problem 
will  be  z?  =  0,  and  x  =  x° •   Now  the  original  problem  can  be 
solved  using  x°  as  the  initial  feasible  point. 
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III.   THE  COMPUTER  PROGRAM  FEASBL 

This  section  will  describe  the  computer  programming  aspects 
of  the  code,  a  complete  set  of  user  instructions  is  provided 
in  section  IV.   Several  programming  languages  were  considered 
for  the  program  including  IBM  360/FORTRAN  IV,  SIMSCRIPT  II. 5, 
and  PL/1.   FORTRAN  IV  was  chosen  because  its  virtually  uni- 
versal availability"  and  widespread  familiarity  will  allow 
easy  installation  at  other  computer  centers,  and  facilitate 
user  modifications  in  the  future. 

FEASBL  is  structured  as  a  straight  line  program,  consist- 
ing of  several  sections;  input,  direction  finding  problem 
building  and  solving,  step  length  determination,  and  output. 
In  deciding  the  size  of  the  problem  that  the  program  would 
be  able  to  handle,  the  size  of  the  direction  finding  linear 
programming  problem  was  the  controlling  factor.   To  insure 
that  round-off  errors  would  not  degrade  the  solution,  the 
inverse  of  the  basis  matrix  is  maintained  in  double  precision. 
FEASBL  is  capable  of  handling  fifty  variables,  with  fifty 
non-linear  and  fifty  linear  constraints  thus  requiring  a 
basis  matrix  that  is  152  by  152,  which  requires  over  184K  of 
core.   Desiring  to  keep  the  total  size  below  300K  if  pos- 
sible, substantial  effort  was  made  to  minimize  the  storage 
requirements  throughout  the  program  both  by  not  storing  any 
value  that  could  be  computed  when  needed,  and  using  the  same 
storage  space  for  multiple  purposes.   The  code  now  requires 
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about  290K  of  core  space  in  addition  to  that  required  for 
user  provided  functions. 

The  information  required  by  the  code  includes  the  problem 
description,  the  initial  starting  point,  and  some  program 
control  variables.   The  problem  description  is  provided  by 
specifying  the  objective  function,  non-linear  and  linear 
constraints,  upper  and  lower  bounds  on  the  variables,  and  the 
requirements  vectors.   In  addition,  the  first  partial  deriva- 
tives of  the  objective  function  and  the  non-linear  constraint 
functions  are  required.   A  prototype  code  was  developed  which 
did  not  require  derivatives,  but  which  severely  limited  the 
form  of  the  functions  that  could  be  used.   It  was  felt,  how- 
ever, that  the  greatly  increased  power  and  flexibility  gained 
by  not  restricting  the  form  of  the  functions  more  than  offset 
the  requirement  to  provide  the  derivatives.   These  functions 
are  provided  in  the  form  of  Fortran  FUNCTION  routines,  thus 
FEASBL  can  handle  any  function  the  user  can  express  as  a 
Fortran  routine.   The  remaining  data  is  input  in  the  Fortran 
NAMELIST  format.   This  method  was  chosen  because  the  unstruc- 
tured format  minimizes  keypunch  errors. 

Each  iteration  begins  by  building  the  direction  finding 
problem.   The  sets  I„(x,e),  IT(x,e),  J  (x,e)  and  J  (x,e)  are 
determined,  and  then  the  direction  finding  linear  programming 
problem  is  built.   Rather  then  store  the  entire  A  matrix, 
only  those  elements  which  depend  on  the  gradients  of  the  ob- 
jective function  and  non-linear  constraints  in  I„(x,e)  referred 
to  as  the  core  are  stored.   The  linear  programming  problem  is 
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solved  using  the  revised  simplex  algorithm  [3],  and  the  basis 
inverse  is  maintained  through  the  product  form.   At  each  step, 
every  non-basic  variable  is  priced,  and  that  variable  with 
the  most  negative  z.  -  c.  is  selected  for  entry.   In  the 
event  of  ties  the  left-most  variable  is  picked.   The  column 
vector  a.,  for  the  entering  variable  is  then  generated,  and 
used  to  compute  v_,  =  B    •  a,  .   The  departing  variable  is  j 
such  that  x,  /y.  ,  is  minimum  for  all  y.  v  >  0.   If  there 

D-j   J  > K  J)K 

are  no  negative  z-  -  c's,  then  the  optimal  solution  has 
been  found  for  the  direction  finding  problem.   After  an  op- 
timal solution  has  been  found,  the  basis  matrix  is  created 
and  inverted  directly  using  the  IBM  subroutine  DMINV,  and 
the  solution  is  checked  for  optimality.   If  the  solution  is 
not  optimal,  then  the  linear  program  is  continued  until  an 
optimal  solution  is  reached.   When  the  optimal  solution  is 
achieved,  the  direction  is  available,  and  the  program  pro- 
ceeds to  find  the  best  step  length.   This  is  accomplished  by 

k     k 
first  finding  Ap,  the  largest  A  such  that  x   +  As_  is  fea- 

K    k 
sible,  and  then  determining  A,  <  Ap  that  maximizes  F(x  +A,s_  ). 

The  quantity  Ap  is  found  by  starting  with  A  =  1  and  increas- 
ing A  by  a  factor  of  2  until  a  AMF  is  found  such  that 

k        k 
x   +  A.IT,  s   is  not  feasible.   The  interval  between  A,IT,  and 
—     NF  —  NF 

the  largest  known  feasible  step  XI.    is  reduced  by  checking 
A  =  ^(^vr  +  AA)  for  feasibility  until  the  interval  is  less 

then  the  minimum  step  length  (STPMIN) .   However,  if  a  feasi- 

k     k       k 
ble  step  length  A  is  found  such  that  F  (x   +  As_  )  <  F  (x  ), 

then  it  is  assumed  that  the  maximum  point  along  the  direction 
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lies  between  these  two  points.   Once  A   has  been  determined 

A,  is  computed  by  the  method  of  regula   falsi.      Since  F  (x)  is 

k  k 

concave  and  is  increasing  at  x   along  the  direction  s_  the 

T  k 

directional  derivative  (V  F(x)   s_)  is  positive  at  x   and  will 

.A. 

be  zero  at  the  maximum  along  s_  .   If  the  directional  deriva- 

k      k 
tive  is  negative  at  x   +  Aps_  ,  then  the  maximum  lies  between 

the  two  points.   A  new  A  is  then  selected  such  that  A  = 
(bX-  -  aA9)  /  (b  -  a)  where:  a  =  V  F(xk  +  A  sk) ;  b  =  V  F(xk 
+  A?s_  );  and  A-  and  A?  bracket  the  optimal  stop  length.   If 
the  directional  derivative  is  zero  at  the  new  point  the  maxi- 
mum has  been  found,  otherwise  we  have  a  new  interval  with 
positive  and  negative  directional  derivatives  at  the  end- 
points,  and  another  A  is  determined.   This  continues  until 
either  the  maximum  is  found,  or  the  interval  becomes  less 

then  the  minimum  specified.   If  the  directional  derivative 

k      k        k      k       k 
is  positive  at  x   +  A  s_  and  F(x   +  A  s_  )  >  F(x  ),  then  the 

maximum  along  s_  is  outside  the  feasible  region,  and  A,  is 

k      k       k 
set  equal  to  Ap.   However,  if  F  (x   +  Aps_  )  <  F  (x  ),  then  the 

function  is  not  concave,  so  a  step  A.,  between  0  and  A„  is 

found  such  that  the  directional  derivative  is  negative,  and 

k      k      k 
it  is  assumed  that  F(x)  is  concave  between  x   and  x   +  A..s_  , 

and  the  maximum  is  sought  as  before. 

Once  a  usable  feasible  direction,  and  step  length  have 

k+ 1    k      k 
been  determined  set  x   '  =  x   +  A,  s_  and  repeat  the  process 

starting  at  x    .   This  process  is  continued  until  either 

no  usable  feasible  direction  can  be  found  or  the  problem 

becomes  unbounded.   If  the  step  length  or  the  percent 
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improvement  in  the  objective  function  become  less  then  user 
supplied  limits  the  program  will  terminate.   The  program  will 
also  stop  after  a  user  specified  number  of  iterations. 
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IV.   USING  FEASBL 

This  section  will  describe  how  to  use  FEASBL,  including 
required  input,  possible  errors,  and  the  normal  output.  An 
example  is  given  illustrating  different  features  of  the  code. 

A.   INPUT 

The  objective  function,  non-linear  constraint  functions, 
and  their  first  partial  derivatives  are  input  as  Fortran 
FUNCTION  routines.   The  linear  constraint  matrix,  require- 
ments vectors,  upper  and  lower  bounds  on  x.,  the  initial  fea- 
sible point,  and  any  other  input  is  read  in  as  data  at 
execution. 

1.   The  FUNCTION  Routines 

Four  Fortran  FUNCTION  routines  are  required,  these 

routines  provide  the  objective  function,  its  first  partial 

derivatives,  the  non-linear  constraints,  and  their  first 

partial  derivatives.   Even  if  the  problem  has  no  non-linear 

constraints,  all  four  FUNCTION  routines  must  be  provided. 

The  four  routines  are: 

FUNCTION  Definition 

FX(  X  )  The  objective  function. 

DELFX(  NX,  X  )      The  partial  derivatives  of 

the  objective  function. 

GIX(  NGI,  X  )      The  non-linear  constraint 

functions . 

DELGIX(NX,NGI,X)    The  partial  derivatives  of 

the  non-linear  constraints, 
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The  arguments  of  the  FUNCTIONS  statements  are 


Argument 
X 
NX 

NGI 


Definition 

The  vector  x. 

Subscript  of  the  partial 
derivative  desired,  i.e. 
9F(x)/9xNX. 

The  non-linear  constraint 
number. 


Each  routine  must  include  a  statement  dimensioning  X(e.g.; 
REALM  X(8),  or  DIMENSION  X(15)  ),  at  least  one  assignment 
statement  with  the  name  of  the  FUNCTION  on  the  left  of  the 
equal  sign  (e.g.;  FX  =  X(l),  or  DELGIX  =  4.0  *  X(2)**3  ), 
and  at  least  one  RETURN  statement.   No  variation  in  the  names 
of  the  FUNCTIONS,  nor  in  the  order  of  the  arguments  is  al- 
lowed.  The  following  example  illustrates  the  use  of  the 
FUNCTION  routines. 


Max: 
Subj  to 


F(x)  =  6x1  +  8x?    x^ 


2     1     2 


a    2  2 

4xl  +  X2 


<  16 


non-linear  constraints 


(4  -  x1)2  +  x22    <    16 


3x1  +  5x~ 


2x1  +  12x? 


0  <  x1    <    1.85 


2.0  <  x. 


<  15 

<  30 


\   linear  constraints 
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For  this  problem  the  following  FUNCTION  routines  are  required: 

a.  The  objective  function 

FUNCTION  FX(  X  ) 

REALM  X(2) 

FX  =  6.0  *  X(l)  +  8.0  *  X(2)  -  X(l)**2  -  X(2)**2 

RETURN 

END 

b.  The  partial  derivatives  of  the  objective  function 

FUNCTION  DELFX(  NX,  X  ) 

REALM  X(2) 

GO  TO  (  101,  102  ) ,  NX 

101  DELFX  =  6.0  -  2.0  *  X(l) 
RETURN 

102  DELFX  =  8.0  -  2.0  *  X(2) 
RETURN 

END 

c.  The  non-linear  constraint  functions 

FUNCTION  GIX(  NGI ,  X  ) 

REALM  X(2) 

GO  TO  (  100,  200  ) ,  NGI 
100  GIX  =  4.0  *  X(l)**2  +  X(2)**2 

RETURN 
200  GIX  =  (  4.0  -  X(l)  )**2  +  X(2)**2 

RETURN 

END 

d.  The  partial  derivatives  of  the  non-linear  con- 
straints 

FUNCTION  DELGIX(  NX,  NGI,  X  ) 

REALM  X(2) 

GO  TO  (  100,  200  ) ,  NGI 

100  GO  TO  (  101,  102  ) ,  NX 

101  DELGIX  =  8.0  *  X(l) 
RETURN 

102  DELGIX  =  2.0  *  X(2) 
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RETURN 

200  GO  TO  (  201,  202  ) ,  NX 

201  DELGIX  =  2.0  *  X(l)  -  8.0 
RETURN 

202  DELGIX  =  2.0  *  X(2) 
RETURN 

END 

The  significant  feature  to  note  is  that  when  a  FUNCTION  is 
called  by  FEASBL  specifying  a  point  X,  and  when  necessary 
the  partial  derivative  desired  NX,  and  the  constraint  number 
NGI ,  the  appropriate  value  is  returned. 
2 .   Data  Input 

The  remainder  of  the  data  needed  is  read  in  from 
data  cards  at  execution  in  NAMELIST  format.   This  format  was 
chosen  to  minimize  keypunch  errors,  and  provide  flexibility 
for  the  user.   The  following  variables  are  required: 

Name  Description 

IXMAX  Number  of  variables. 

IGCNT  Number  of  non-linear  constraints 

ILCNT  Number  of  linear  constraints. 

The  following  variables  are  optional,  and  if  not  provided  by 
the  user  will  take  on  the  default  values  listed. 

Name  Default  Description 

GL(I,J)      0.0  Linear  constraint  matrix  element 

XMIN(I)      0.0  Lower  bound  on  x, 

XMAX(I)      =1075  Upper  bound  on  x. 

BC(I)        0.0      The  b.  for  the  non-linear  con- 
straints . 
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Name 

Default 

BL(I) 

0.0 

X(I) 

0.0 

EPSLON 

0.1 

STPMIN 

10-3 

PERCNT 

10"2 

UMAX 
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Description 

The  b.  for  the  linear  constraints, 

The  initial  feasible  point. 

The  initial  anti- jamming  constant, 

The  minimum  step  length. 

The  minimum  objective  function 
improvement . 

The  maximum  number  of  iterations. 


The  first  data  card  begins  with  an  ampersand  (§)  in  column  two 
followed  by  the  word  INPUT  with  no  space  (  §INPUT"  )  .   The 
values  are  then  listed  as  "variable  name  =  FORTRAN  constant" 
i.e.  IXMAX  =  2,  GL(1,1)  =  3.0,  the  last  data  item  is  followed 
by  an  ampersand  and  the  word  "END"  (  §END  ) .   If  more  then 
one  card  is  required,  each  card  must  begin  in  column  2  or 
more,  and  the  §END  will  appear  on  the  last  card  only.   For 
the  above  example  the  data  cards  would  look  like: 

SINPUT  IXMAX  =  2,  IGCNT=2,  ILCNT=2,  GL(1,1)=3.0,  GL(1,2)=5.0, 
GL(2,1)=2.0,  GL(2,2)=12.0,  XMAX (1) = 1 . 85 ,  XMIN(2)=2.0, 
BC(1)=16.0,  BC(2)=16.0,  BL(1)=15.0,  BL(2)=30.0,  ITMAX=10, 
EPSLON=0.01,  STPMIN=1.0E-04,  X(l)=1.0,  X(2)=2.0,  SEND 
Here  the  order  is  unimportant,  and  all  the  variables  need  not 
be  listed.   In  this  example  the  default  values  for  XMIN(l), 
XMAX (2),  and  PERCNT  will  be  used  since  they  are  not  listed 
on  the  data  cards. 
3.   Data  Limits 

Input  errors  will  result  if  the  variable  inputs  are 
outside  the  following  limits: 
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1  <  IXMAX  <  50 
0  <  IGCNT  <  50 
0  <  ILCNT  <  50 

XMIN(I)  <  XMAX(I) 

To  prevent  stopping  prematurely,  EPSLON  should  be  significantly 
larger  than  STPMIN. 

After  all  the  data  is  input,  the  initial  point  is 
checked  for  feasibility.   If  the  point  is  not  feasible  the 
program  is  terminated. 

4.  The  Program  CHECK 

FEASBL  requires  290K  of  core  space  in  addition  to 
that  required  for  the  FUNCTION  routines.   Most  problems  will 
require  less  then  320K  and  less  then  4  minutes  to  run.   As 
an  aid  for  checking  the  input  data  and  FUNCTION  routines,  the 
program  CHECK  can  be  used.   CHECK  requires  exactly  the  same 
input  as  FEASBL,  and  performs  the  same  input  error  checks. 
In  addition  CHECK  will  determine  the  feasibility  of  the  ini- 
tial point  and  print  out  the  value  of  the  objective  function, 
all  the  constraints,  and  the  value  of  all  the  partial  deriva- 
tives of  the  objective  function  and  non-linear  constraints. 
By  using  CHECK  the  user  can  quickly  verify  the  input,  since 
CHECK  can  be  run  normally  as  a  class  B  or  C  job,  before 
using  FEASBL,  which  is  a  class  J  job. 

5 .  System  Cards 

The  JCL  cards  required  to  use  FEASBL  or  CHECK  depend 
on  the  way  that  the  FUNCTION  routines  are  input.   The  routines 
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can  be  input' as  either  Fortran  source  decks  and  compiled  as 
part  of  the  job,  or  the  routines  can  be  precompiled  and  in- 
put as  object  decks.   An  example  of  each  method  follows. 

a.   FEASBL  with  Source  Decks 

//    standard  OS  JOB  card       ,TIME= (03 , 59) 
//  EXEC  FORTCLG, REGION. GO=300K 
//FORT.SYSIN  DD  * 

(FUNCTION  source  decks) 

/* 

//LINK.NLPLIB  DD  DSNAME=F0600 . NLPLIB ,UNIT=2321 , 
//    VOL=SER=CEL001,DISP=SHR,LABEL=(, , ,IN) 
//LINK.SYSIN  DD  * 

INCLUDE  NLPLIB (FEASBL, DARRAY) 

/* 

//GO.SYSIN  DD  * 
SINPUT  IXMAX-2,  IGCNT=2,  ...  ,  SEND 
/* 
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b.   FEASBL  with  Object  Decks 

//    standard  OS  JOB  card      ,TIME= (03 , 59) 
//  EXEC  LGO, REGION. GO=300K 

//LINK.NLPLIB  DD  DSNAME=F0600 . NLPLIB ,UNIT=2321 , 
//    VOL=SER=CEL001,DISP-SHR,LABEL=C, , ,IN) 
//LINK.SYSIN  DD  * 

INCLUDE  NLPLIB (FEASBL, DARRAY) 
(FUNCTION  object  decks) 

/* 

//GO.SYSIN  DD  * 

SINPUT  IXMAX=2,  IGCNT=2,  ...  ,  X(2)=2.0,  SEND 
/* 


30 


c.   CHECK  with  Source  Decks 

//    standard  OS  JOB  card 

//  EXEC  FORTCLG, REGION. GO=100K 

//FORT.SYSIN  DD  * 

(FUNCTION  source  decks) 

/* 

//LINK.NLPLIB  DD  DSNAME=F0600 . NLPLIB ,UNIT=2321 , 
//   VOL=SER=CEL001,DISP=SHR,LABEL=(, , ,IN) 
//LINK.SYSIN  DD  * 

INCLUDE  NLPLIB (CHECK) 

/* 

//GO.SYSIN  DD  * 

§INPUT  IXMAX=2,  IGCNT=2,  ...  ,X(2)=2.0,  SEND 
/* 
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d.   CHECK  with  Object  Decks 

//      standard  OS  JOB  card 

//  EXEC  LGO, REGION. GO=100K 

//LINK.NLPLIB  DD  DSNAME=F0600 . NLPLIB,UNIT=2321 , 

//   VOL=SER=CEL001,DISP=SHR,LABEL=(, , ,IN) 

//LINK.SYSIN  DD  * 

INCLUDE  NLPLIB(CHECK) 

(FUNCTION  objective  decks) 

/* 

//GO.SYSIN  DD  * 

SINPUT  IXMAX=2,  IGCNT=2,  ...  ,  X(2)=2.0,  SEND 
/* 

A  set  o£  object  decks  can  be  produced  either  after 
checking  the  FUNCTION  routines  by  executing  a  FORTCD  job,  or 
in  conjunction  with  the  program  CHECK  by  replacing  //  EXEC 
FORTCLG, REGION. GO=100K  with  //  EXEC  FORTCLGP ,PARM. F0RT= 
'DECK* , REGION. GO=100K. 

B.   SAMPLE  PROGRAM 

The  FUNCTION  routines,  input  cards,  and  system  cards 
required  to  use  FEASBL  with  source  decks  are  shown  for  the 
problem: 


Max 


2     2 
F(x)  =  6x1  +  8x2  -  x?  -  x» 


Subj  to:     4x1  +  x-       <    16 
(4  -  xx)2  +  x22    <    16 
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3x1  +  5x„      <  15 

2x1  +  12x2      <  30 

0 .  <  x-  <  1.85 

2.0  <  x   <  °° 

A  sample  o£  the  output  o£  both  CHECK  and  FEASBL  for  this 
problem  is  provided.   The  next  two  pages  are  the  cards  re- 
quired for  executing  FEASBL  with  the  sample  problem.   This 
is  followed  by  the  output  of  CHECK  and  then  FEASBL  for  the 
sample  problem. 
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//SAMPLE01    JOB    (1115, 0601, ROY2) ,' SAMPLE    JOB •  , T I  ME  =  ( 03 , 59 ) 
//    EXEC    FORTCLGi REGION. G0=300K 
//FOF-T.SYSIN    DD    * 

FUNCTION    FX(    X    ) 
C 
C  THIS    FUNCTION    PROVIDES    THE    OBJECTIVE    FUNCTION. 


C 

c 
c 
c 


REAL*4    X(2) 
C 

C  THIS    IS    THE    OBJECTIVE    FUNCTION 


FX    =    6.0    *    X(l)     +    8.0    *    X12)    -    X(l)**2    -    X(2)**2 

RETURN 
END 
C 

C 

C  ***************************************** 

c 
c 

FUNCTION    DELFX(     NX,    X     ) 
C 

C  THIS    FUNCTION    PROVIDES    THE    FIRST    PARTIAL    DERIVATIVES 

C  OF    THE    OBJECTIVE    FUNCTION. 

C 

REAL*4    X(2) 

GO    TO    (     101,     102     ) ,    NX 
C 

C  THIS    IS    THE    PARTIAL    WITH    RESPECT    TO    X(l). 

C 

101  DELFX    =    6.0   -    2.0    *    X(l) 
RETURN 

C 

C  THIS    IS    THE    PARTIAL    WITH    RESPECT    T0    X(2) 

C 

102  DELFX    =    8.0    -    2.0    *    X(2) 
RETURN 

END 
C 

c 

f~  3p  ^f  ^  ','  i»  »,»  5  '  f  *  ^C  i|«  Ip  •$*  *v  -P  ***  X*  "V  *  vr  f*  *f*  **  **f*  *  *t*  t  3jS  3p  SJs  3|6  3jK  Sp  3p  3p  ^*  3J*  ^6  3j*  *F  «»* 

C 

C 

FUNCTION    GIX(    NGI,     X    ) 
C 

C  THIS    FUNCTION    PROVIDES    THE    NON-LINEAR    CONSTRAINTS. 

C 

REAL*4    X(2) 
C 

C  GO    TO    THE    PROPER    CONSTRAINT. 

C 

GO    TO    (     100,    200     ) ,    NGI 
C 

C  THIS     IS    NON-LINEAR    CONSTRAINT    NUMBER    1. 

C 

100    GIX    =    4.0    *    X(l)**2    +    X(2}**2 

RETURN 
C 

C  THIS    IS    NON-LINEAR    CONSTRAINT    NUMBER    2. 

C 

200    GIX    =     (    4.0    -    X{1)     )**2    +    X(2)**2 

RETURN 

END 
C 
C 

Q  ****  *******  *****  *****************  ******** 

c 

C 

FUNCTION    DELGIX(     NX,    NGI,     X    ) 
C 

C  THIS    FUNCTION    PROVIDES    THE    FI^ST    PARTIAL    DERIVATIVES 

C  OF    THE    NON-LINEAR    CONSTRAINT    FUNCTIONS. 

C 
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REAL*4  X(2) 
C 

C  GO  TO  THE  CONSTRAINT  SPECIFIED. 

C 

GO  TO  (  100,  200  ),NGI 
C 

C      THIS  IS  CONSTRAINT  NUMBER  ONE,  WHICH  PARTIAL  IS  WANTED 
C 

100  GO  TO  (  101,  102  ),  NX 
C 

C  THIS  IS  THE  PARTIAL  WITH  RESPECT  TO  X(l). 

C 

101  DELGIX  =  8.0  *  X(  1) 
RETURN 

C 

C  THIS  IS  THE  PARTIAL  WITH  RESPECT  TO  X(2) 

C 

102  DELGIX    =    2.0    *    X(2) 
RETURN 

C 

C      THIS  IS  CONSTRAINT  NUMBER  TWO,  WHICH  PARTIAL  IS  WANTED 

C 

200  GO  TO  (  201,  202  ),  NX 
C 

C  THIS  IS  THE  PARTIAL  WITH  RESPECT  TO  X(l) 

C 

201  DELGIX  =  2.0  *  X(l)  -  6.0 
RETURN 

C 

C      THIS  IS  THE  PARTIAL  WITH  RESPECT  TO  X(2) 

C 

202  DELGIX  =  2.0  *  X(2) 
RETURN 

END 

//LINK.NLPLIB  DD  DSN=F0600. NLPL I B ,UN I T=232 1 , 

//        VCL=SER=CEL001 , DI S P=SHR, LA BEL= ( , , , IN ) 

//LINK.SYSIN  DD  * 

INCLUDE  NLPLIB(FEASBL,DARRAY) 

//GO.SYSIN  DD  * 
&INPUT  IXMAX=2,  IGCNT=2,  ILCNT=2,  GL(1,1)=3.0,  GL(1,2)=5.0, 
GL(2,1)=2.0,  GL(2,2)  =  12.C,  XMAX( 1 1  =  1  .  85 ,  XMIN(2)=2.0, 
BC(1)=16.0,  3C(2)=16.0,  BL(1)=15.0,  BL(2)=30.0,  ITMAX=10, 
EPSLCN=0.01,  STPMIN=1.0E-04f  X(l)=1.0,  X(2)=2.0,  &END 
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SAMPLE      OUTPUT     FROM     CHECK 


NON-LINEAR  PROGRAMMING  PROBLEM 


INITIAL  SET  UP: 

NUMBER  OF  VARIABLES:  2 

NUMBER  OF  NON-LINEAR  CONSTRAINTS:     2 

NUMBER  OF  LINEAR  CONSTRAINTS:  2 

THE  LINEAR  CONSTRAINT  MATRIX  IS: 

(  li  1);     3.00030  (  1,  2);     5.00000 


(  2.  1); 


2.00000 


(  2.  2); 


12.0000 


THE  BOUNDS  ON  THE  VARIABLES  ARE: 

0.0  <  XI  1)  <      1.85000     : 

NCN-LINEAR  CONSTRAINTS  REQUIREMENTS  VECTOR,  BCIII: 

(  l);         I6.0000    (  2):         I6.0000 

LINEAR  CONSTRAINTS  REQUIREMENTS  VECTOR,  BL ( I ) : 

t  II;  15.0000     (  2):  30.0000 

THE  INITIAL  FEASIBLE  STARTING  POINT  IS: 


X(  1); 

ANTI  JAMMING  CONSTANT 
MINIMUM  STEP  LENGTH 
MINIMUM  IMPROVEMENT 
MAXIMUM   «  ITERATIONS 


1.00000 


,xi    2): 

O.IOOOOOE-Oi 

0.  100000E-03 

O.IOOOOOE-OI 

10 


INPUT    COMPLETE.        INITIAL    VALUE    OF    OBJECTIVE    FUNCUON: 


2.00000 


<    X(    2)    < 


0.723701E    76 


2.00000 


,X( 


1.700000E    01 


VALUE    OF    NON-LINEAR    CONSTRAINT    NUMBER,       1,     IS:        8.0000000E    00 

VALUE    OF    NON-LINEAR    CONSTRAINT     NUMBER,       2,     IS:        1.3000000E    01 

THE    VALUE    OF       LINEAR    CONSTRAINT    NUMBER        1,     IS:  1.300000E    01 

THE    VALUE    OF       LINEAR    CONSTRAINT    NUMBER       2,      IS:  2.600000E    01 

THE    BCUNOS,     ANO    ALL    THE    NCN-LINEAR    AND    LINEAR    CONSTRAINTS    HAVE    BEEN    CHECKED. 

NEXT    THE    VALUES    OF    THE    PARTIAL    DERIVATIVES    OF    THE    OBJECTIVE    FUNCTION    AND    THE    CONSTRAINTS    ARE    LISTED. 

THE    VALUE    OF    OELFXt     1,     XI     IS:  <,.0O0000E    00 

THE    VALUE    OF    DELFXI     2,     X)      IS:  4.000000E     00 

THE    VALUE    OF    DELGIXI     1,     1,    X)     IS:  8.000000E    00 

THE    VALUE    OF    DELGIXI     2,     1,    XI     IS:  <..000000E    00 

THE    VALUE    OF    DELGIXI     1,     2,     X)     IS:       -6.000000E     00 

THE    VALUE    OF    DELGIXI     2,     2,     X)     IS:  <,.000000E     00 
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SAMPLE      OUTPUT     FROM     FEASBL 


NON-LINEAR  PROGRAMMING  PROBLEM 


INITIAL  SET  UP: 

NUMBER  OF  VARIA3LES:  2 

NUMBER  OF  NON-LINEAR  CONSTRAINTS:     2 

NUMBER  OF  LINEAR  CONSTRAINTS:         2 

THE  LINEAR  CONSTRAINT  MATRIX  IS: 

(  1,  1);     3.00000  (  1,  2);     5.00000 


(  2,  1); 


2.00000 


(  2,  2); 


12.0000 


THE  BOUNDS  ON  THE  VARIABLES  ARE: 

0.0  <  X(  1)  <      1.85000     ;  2.00000      <  XI  2)  < 

NON-LINEAF  CONSTRAINTS  REQUIREMENTS  VECTOR.  BC(I|: 
(  1):  16.0000     (  2):  16.0000     I 

LINEAR  CCNSTRAINTS  REQUIREMENTS  VECTOR,  6L I  I  I  : 
.(  1);  15.0000     (  2);  30.0000     I 

THE  INITIAL  FEASIBLE  STARTING  POINT  IS: 

X(  1);  1.00000    ,X(  21;  2.00000     ,X( 

ANT!  JAMMING  CONSTANT:  0.100000E-01 

MINIMUM   STEP   LENGTH:  0.100000E-03 

MINIMUM    IMPROVEMENT:  0.100000E-01 

MAXIMUM  *    1TERATICNS:  10 

INPUT  COMPLETE.   INITIAL  VALUE  OF  OBJECTIVE  FUNCTION:     1.700000E  01 


0.723701E    76 


,TERF(X^N=NUyBl!8S7500E    01S        IMPROVEMENT    THIS     ITERATION: 


PROBLEM    SOLUTION 
11.0294    PERCENT. 


THE    POINT    X{     )     IS: 

X(     11,  1.250000 


XI     21  , 


2.250000 


XI 


ITERATION    NUMBER:  2 

FIX)     a  1.922217E    01;       IMPROVEMENT    THIS     ITERATION: 


1.8393    PERCENT. 


THE    POINT    X(     I     IS: 

X(  1),     1.666625     ;  X(  21,     2.000024     ;  X( 
OPTIMAL  SOLUTION: 
THE  OPTIMAL  VALUE  OF  THE  OBJECTIVE  FUNCTION  IS:    1.922217E  01 

THE  PCINT  XI  1  IS: 

XI  II,      1.666625      ;  X(  2),     2.000024      ;  X( 
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V.   EVALUATION 

To  validate  the  code,  and  evaluate  its  performance, 
several  non-linear  problems  were  solved  with  FEASBL.   These 
problems  were  selected  to  represent  the  range  of  problems 
FEASBL  was  designed  to  solve.   The  problems  included  uncon- 
strained optimization  of  a  concave  objective  function,  con- 
strained optimization  of  concave  objective  functions  some 
with  the  unconstrained  maxima  interior  and  some  exterior 
to  the  feasible  regions.   An  unbounded  problem  was  also 
solved. 

In  each  case  the  optimum  solution  was  obtained  in  less 
then  five  seconds  of  cpu  time.   These  problems  were  all  small 
in  order  to  facilitate  manual  verification  of  the  solutions. 
Larger  problems  will  take  correspondingly  more  time,  but  the 
code  appears  to  be  quite  efficient.   In  most  problems  the 
optimal  solution  was  obtained  after  two  iterations.   In  one 
case,  the  solution  took  seven  iterations. 

It  was  found  that  the  value  of  STPMIN  is  critical  to  the 

performance  of  FEASBL.   If  STPMIN  is  set  too  small,  when  the 

interval  between  X„  and  X..„  is  calculated,  the  round-off 

F      NF  ' 

error  may  exceed  the  value  of  STPMIN.   This  results  in  the 
code  looping  indefinitely  trying  to  reduce  the  interval  beyond 
machine  accuracy  leading  to  a  program  completion  due  to  ex- 
cess time.   If  this  should  occur  the  problem  should  be  re- 
started at  the  last  point  with  a  larger  value  for  STPMIN. 
With  step  lengths  on  the  order  of  1.0  the  interval  cannot  be 
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determined  to  an  accuracy  greater  then  10   .   With  larger 
step  lengths  this  value  will  increase.   An  initial  value 
of  10    or  10"3  is  recommended  for  STPMIN.   This  value 
can  be  reduced  in  subsequent  runs  of  FEASBL  if  a  more  pre- 
cise solution  is  desired. 

The  initial  value  of  EPSLON  is  not  as  critical  as  STPMIN. 
EPSLON  should  be  greater  than  STPMIN  but  generally  not  greater 

than  10   .   Experience  with  FEASBL  indicates  that  values  on 

- 1      -  2 
the  order  of  10    or  10   perform  well.   If  a  large  value  is 

used,  on  the  order  of  1.0  or  greater,  and  an  optimal  point 

is  found,  this  should  be  verified  using  a  smaller  value  for 

EPSLON,  as  the  code  may  stop  early  with  large  values  of 

EPSLON. 

Performance  was  uneffected  by  variations  in  PERCNT.   Large 
values  may  cause  the  program  to  terminate  sooner.   Values  less 
then  10    are  equivalent  to  zero  because  of  inherent  machine 
accuracy . 

An  attempt  was  made  to  solve  a  problem  in  which  all  the 
constraints  were  equality  constraints,  thus  the  feasible 
region  had  no  strictly  interior  points.   The  program  would 
terminate  immediately  without  any  significant  movement.  This 
phenomenon  was  anticipated  since  the  underlying  theory  assumes 
that  the  feasible  region  contains  at  least  one  strictly  in- 
terior point  to  the  set  of  non-linear  constraints  [1]. 
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